{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "name": "python3",
   "language": "python"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.14",
   "mimetype": "text/x-python",
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "pygments_lexer": "ipython3",
   "nbconvert_exporter": "python",
   "file_extension": ".py"
  },
  "kaggle": {
   "accelerator": "nvidiaTeslaT4",
   "dataSources": [],
   "dockerImageVersionId": 30776,
   "isInternetEnabled": true,
   "language": "python",
   "sourceType": "notebook",
   "isGpuEnabled": true
  },
  "colab": {
   "provenance": [],
   "gpuType": "T4"
  },
  "accelerator": "GPU"
 },
 "nbformat_minor": 5,
 "nbformat": 4,
 "cells": [
  {
   "cell_type": "code",
   "source": [
    "import torch\n",
    "from collections import OrderedDict\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io\n",
    "#from pyDOE import lhs#拉丁超立方抽样\n",
    "from torch import nn\n",
    "import time\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)  # 只忽略UserWarning类型的警告"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:44.161933Z",
     "iopub.execute_input": "2024-09-24T08:23:44.162356Z",
     "iopub.status.idle": "2024-09-24T08:23:45.981338Z",
     "shell.execute_reply.started": "2024-09-24T08:23:44.162314Z",
     "shell.execute_reply": "2024-09-24T08:23:45.980437Z"
    },
    "trusted": true,
    "id": "af4b048c0e49992d",
    "ExecuteTime": {
     "end_time": "2024-12-03T13:31:12.860520Z",
     "start_time": "2024-12-03T13:31:11.522127Z"
    }
   },
   "id": "af4b048c0e49992d",
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "if torch.cuda.is_available():\n",
    "    device = torch.device('cuda')\n",
    "else:\n",
    "    device = torch.device('cpu')\n",
    "def setup_seed(seed):\n",
    "    torch.manual_seed(seed)\n",
    "    torch.cuda.manual_seed_all(seed)\n",
    "    np.random.seed(seed)\n",
    "    #  random.seed(seed)\n",
    "    torch.backends.cudnn.deterministic = True\n",
    "\n",
    "setup_seed(124)\n",
    "dtype = torch.double"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:45.986195Z",
     "iopub.execute_input": "2024-09-24T08:23:45.988477Z",
     "iopub.status.idle": "2024-09-24T08:23:46.049102Z",
     "shell.execute_reply.started": "2024-09-24T08:23:45.988430Z",
     "shell.execute_reply": "2024-09-24T08:23:46.048271Z"
    },
    "trusted": true,
    "id": "9e2f9d35aa33fbb9",
    "ExecuteTime": {
     "end_time": "2024-12-03T13:31:12.876127Z",
     "start_time": "2024-12-03T13:31:12.861537Z"
    }
   },
   "id": "9e2f9d35aa33fbb9",
   "execution_count": 2,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#求导\n",
    "def gradients(outputs, inputs):\n",
    "    return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True)\n",
    "#类型转换\n",
    "def to_numpy(input):\n",
    "    if isinstance(input, torch.Tensor):\n",
    "        return input.detach().cpu().numpy()\n",
    "    elif isinstance(input, np.ndarray):\n",
    "        return input\n"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.052132Z",
     "iopub.execute_input": "2024-09-24T08:23:46.052502Z",
     "iopub.status.idle": "2024-09-24T08:23:46.057747Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.052456Z",
     "shell.execute_reply": "2024-09-24T08:23:46.056804Z"
    },
    "trusted": true,
    "id": "e501d8818e912d51",
    "ExecuteTime": {
     "end_time": "2024-12-03T13:31:12.891687Z",
     "start_time": "2024-12-03T13:31:12.879135Z"
    }
   },
   "id": "e501d8818e912d51",
   "execution_count": 3,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#构建网络layers\n",
    "def Dnnlayers(input_layers,output_layers,hidden_layers,neural):\n",
    "    layers = []\n",
    "    for i in range(hidden_layers + 2):\n",
    "        if i == 0:\n",
    "            layers.append(input_layers)\n",
    "        elif i == hidden_layers + 1:\n",
    "            layers.append(output_layers)\n",
    "        else:\n",
    "            layers.append(neural)\n",
    "    return layers\n",
    "#初值\n",
    "def IC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n",
    "    N = x.shape[0]\n",
    "    rho_init = np.zeros((x.shape[0], 1))\n",
    "    u_init = np.zeros((x.shape[0], 1))\n",
    "    p_init = np.zeros((x.shape[0], 1))\n",
    "\n",
    "    # rho, p - initial condition\n",
    "    for i in range(N):\n",
    "        if (x[i, 1] <= 0.0):\n",
    "            rho_init[i] = crhoL\n",
    "            u_init[i] = cuL\n",
    "            p_init[i] = cpL\n",
    "        else:\n",
    "            rho_init[i] = crhoR\n",
    "            u_init[i] = cuR\n",
    "            p_init[i] = cpR\n",
    "    U_ic = np.hstack([rho_init, u_init, p_init])\n",
    "\n",
    "    return U_ic\n",
    "#外推边界\n",
    "def BC(x, crhoL, cuL, cpL, crhoR, cuR, cpR):\n",
    "    N = x.shape[0]\n",
    "    rho_bc = np.zeros((x.shape[0], 1))\n",
    "    u_bc = np.zeros((x.shape[0], 1))\n",
    "    p_bc = np.zeros((x.shape[0], 1))\n",
    "\n",
    "    # rho, p - initial condition\n",
    "    for i in range(N):\n",
    "        if (x[i, 1] <= 0.0):\n",
    "            rho_bc[i] = crhoL\n",
    "            u_bc[i] = cuL\n",
    "            p_bc[i] = cpL\n",
    "        else:\n",
    "            rho_bc[i] = crhoR\n",
    "            u_bc[i] = cuR\n",
    "            p_bc[i] = cpR\n",
    "    U_bc = np.hstack([rho_bc, u_bc, p_bc])\n",
    "\n",
    "    return U_bc\n",
    "\n",
    "\n",
    "def x_left_right(x_int, dx):\n",
    "    N = x_int.shape[0]\n",
    "    xL = np.zeros((N, 2))\n",
    "    xR = np.zeros((N, 2))\n",
    "\n",
    "    for i in range(N):\n",
    "        xL[i, 1] = x_int[i, 1] - dx\n",
    "        xL[i, 0] = x_int[i, 0]\n",
    "        xR[i, 1] = x_int[i, 1] + dx\n",
    "        xR[i, 0] = x_int[i, 0]\n",
    "\n",
    "    return xL, xR\n"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.060253Z",
     "iopub.execute_input": "2024-09-24T08:23:46.061124Z",
     "iopub.status.idle": "2024-09-24T08:23:46.074353Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.061075Z",
     "shell.execute_reply": "2024-09-24T08:23:46.073511Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:14.962601Z",
     "start_time": "2024-11-13T03:05:14.951807Z"
    },
    "id": "79e6a971aae89fdf"
   },
   "id": "79e6a971aae89fdf",
   "execution_count": 4,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "class DNN(nn.Module):\n",
    "    def __init__(self, layers):\n",
    "        #继承父类\n",
    "        super(DNN, self).__init__()\n",
    "\n",
    "        #depth of network\n",
    "        self.depth = len(layers) - 1\n",
    "        #print(self.depth)\n",
    "\n",
    "        #activation of network\n",
    "        self.activation = nn.Tanh()\n",
    "\n",
    "        #create the neural network\n",
    "        layers_list = list()\n",
    "        for i in range(self.depth - 1):\n",
    "            layers_list.append(\n",
    "                ('layer_%d' % i, nn.Linear(layers[i], layers[i+1]))#create each network\n",
    "            )\n",
    "            layers_list.append(\n",
    "                ('activation_%d' % i, self.activation)  #create each activation of network\n",
    "            )\n",
    "        layers_list.append(\n",
    "            ('layer_%d' % (self.depth - 1), nn.Linear(layers[-2], layers[-1]))\n",
    "            #last network do not have activation\n",
    "        )\n",
    "        #创建一个有序字典，其中包含了从 layers_list 中获得的键值对，这在需要保持元素顺序的场景（如神经网络层的顺序）中非常有用。\n",
    "        layerDict = OrderedDict(layers_list)\n",
    "\n",
    "        #deploy layers\n",
    "        self.layers = nn.Sequential(layerDict)\n",
    "        #print(self.layers)\n",
    "\n",
    "    #forword network,output the result of network\n",
    "    def forword(self, x):\n",
    "        out = self.layers(x)\n",
    "        return out"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.075428Z",
     "iopub.execute_input": "2024-09-24T08:23:46.075700Z",
     "iopub.status.idle": "2024-09-24T08:23:46.086317Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.075671Z",
     "shell.execute_reply": "2024-09-24T08:23:46.085503Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:14.978642Z",
     "start_time": "2024-11-13T03:05:14.964606Z"
    },
    "id": "2b5f4f99b1b8181f"
   },
   "id": "2b5f4f99b1b8181f",
   "execution_count": 5,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "class Pinns():\n",
    "    def __init__(self, layers):\n",
    "        self.dnn = DNN(layers).double().to(device)\n",
    "    #forword network,output the result of network\n",
    "    def forword(self, x):\n",
    "        y = self.dnn.forword(x)\n",
    "        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n",
    "        return rho, u, p\n",
    "    #function loss,残差点训练\n",
    "    def loss_pde(self, x_int):\n",
    "        y = self.dnn.forword(x_int)\n",
    "        rho, u, p = y[:, 0:1], y[:, 1:2], y[:, 2:]\n",
    "        #gamma = self.gamma\n",
    "\n",
    "        U2 = rho * u\n",
    "        U3 = 0.5 * rho * u ** 2 + p / (1.4 - 1.0)#0.4 #gamma - 1 = 0.4\n",
    "        # F1 = U2\n",
    "        F2 = rho * u ** 2 + p\n",
    "        F3 = u * (U3 + p)\n",
    "\n",
    "        # Gradients and partial derivatives\n",
    "        drho_g = gradients(rho, x_int)[0]\n",
    "        rho_t, rho_x = drho_g[:, :1], drho_g[:, 1:]\n",
    "\n",
    "        dU2_g = gradients(U2, x_int)[0]\n",
    "        U2_t, U2_x = dU2_g[:, :1], dU2_g[:, 1:]\n",
    "        dU3_g = gradients(U3, x_int)[0]\n",
    "        U3_t, U3_x = dU3_g[:, :1], dU3_g[:, 1:]\n",
    "        dF2_g = gradients(F2, x_int)[0]\n",
    "        F2_t, F2_x = dF2_g[:, :1], dF2_g[:, 1:]\n",
    "        dF3_g = gradients(F3, x_int)[0]\n",
    "        F3_t, F3_x = dF3_g[:, :1], dF3_g[:, 1:]\n",
    "\n",
    "        f = (((rho_t + U2_x) ) ** 2).mean() + \\\n",
    "            (((U2_t + F2_x) ) ** 2).mean() + \\\n",
    "            (((U3_t + F3_x) ) ** 2).mean()  #\n",
    "\n",
    "        return f\n",
    "    #initial loss function\n",
    "    def loss_ic(self, x_ic, U_ic):\n",
    "        y_ic = self.dnn.forword(x_ic)\n",
    "        rho_ic_nn, u_ic_nn, p_ic_nn = y_ic[:, 0:1], y_ic[:, 1:2], y_ic[:, 2:]\n",
    "        loss_ics = ((rho_ic_nn - U_ic[:,0:1]) ** 2).mean()+\\\n",
    "                   ((u_ic_nn - U_ic[:,1:2]) ** 2).mean()+\\\n",
    "                   ((p_ic_nn - U_ic[:,2:]) ** 2).mean()\n",
    "        return loss_ics\n",
    "    #boundary condition\n",
    "    def loss_bc(self, x_bc, U_bc):\n",
    "        y_bc = self.dnn.forword(x_bc)\n",
    "        rho_bc_nn, u_bc_nn, p_bc_nn = y_bc[:, 0:1], y_bc[:, 1:2], y_bc[:, 2:]\n",
    "        loss_bcs = ((rho_bc_nn - U_bc[:, 0:1]) ** 2).mean() + \\\n",
    "                   ((u_bc_nn - U_bc[:, 1:2]) ** 2).mean() + \\\n",
    "                   ((p_bc_nn - U_bc[:, 2:]) ** 2).mean()\n",
    "        return loss_bcs\n",
    "    def loss_rh(self, x_int, x_int_l, x_int_r):\n",
    "        #all time and all x\n",
    "        y = self.dnn.forword(x_int)\n",
    "        rho, u, p = y[:, 0], y[:, 1], y[:, 2]\n",
    "        e = p/0.40+0.5*rho*u*u\n",
    "        #all time and x+dx\n",
    "        y_r = self.dnn.forword(x_int_r)\n",
    "        rhor, ur, pr = y_r[:, 0], y_r[:, 1], y_r[:, 2]\n",
    "        er = pr/0.40+0.5*rhor*ur*ur\n",
    "        #all time and x-dx\n",
    "        y_l = self.dnn.forword(x_int_l)\n",
    "        rhol, ul, pl = y_l[:, 0], y_l[:, 1], y_l[:, 2]\n",
    "        el = pl/0.40+0.5*rhol*ul*ul\n",
    "\n",
    "        du_g = gradients(u, x_int)[0]\n",
    "        u_t, u_x = du_g[:, :1], du_g[:, 1:]\n",
    "        d = 1 / (0.1 * (abs(u_x) - u_x) + 1)\n",
    "        '''\n",
    "        torch.clamp 函数用于将张量中的元素限制在指定范围内\n",
    "        '''\n",
    "        eta = torch.clamp(abs(p - pl) - 0.1, min=0) * torch.clamp(abs(u - ul) - 0.1, min=0)\n",
    "\n",
    "        loss_rh = (((rhor / rhol - (6 * pr + pl) / (6 * pl + pr)) * (ur - ul) * eta) ** 2).mean() + \\\n",
    "                  ((((ur - ul) ** 2 - 2 / rhor * (pr - pl) ** 2 / (0.4 * pr + 2.4 * pl)) * eta) ** 2).mean()\n",
    "\n",
    "        return loss_rh\n",
    ""
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.087638Z",
     "iopub.execute_input": "2024-09-24T08:23:46.088500Z",
     "iopub.status.idle": "2024-09-24T08:23:46.103784Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.088466Z",
     "shell.execute_reply": "2024-09-24T08:23:46.102979Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.009727Z",
     "start_time": "2024-11-13T03:05:14.981648Z"
    },
    "id": "5a835b6157600a38"
   },
   "id": "5a835b6157600a38",
   "execution_count": 6,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#网格构建\n",
    "num_x = 201\n",
    "num_t = 51\n",
    "Tstart = 0.0\n",
    "Tend = 0.5\n",
    "Xstart = -1.0\n",
    "Xend = 1.0\n",
    "x = np.linspace(Xstart, Xend, num_x)\n",
    "t = np.linspace(Tstart, Tend, num_t)\n",
    "dx = (Xend - Xstart)/(num_x-1)\n",
    "t_grid, x_grid = np.meshgrid(t, x)\n",
    "T = t_grid.flatten()[:, None]\n",
    "X = x_grid.flatten()[:, None]\n",
    "X_star = np.hstack((T.flatten()[:, None], X.flatten()[:, None]))\n",
    "# 时空边界\n",
    "lb = X_star.min(0)  # 时空下界，\n",
    "ub = X_star.max(0)  # 时空上界\n",
    "print('lb is :',lb, ', ub is :', ub)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.413829Z",
     "iopub.execute_input": "2024-09-24T08:23:46.414437Z",
     "iopub.status.idle": "2024-09-24T08:23:46.432210Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.414400Z",
     "shell.execute_reply": "2024-09-24T08:23:46.431284Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.024922Z",
     "start_time": "2024-11-13T03:05:15.011889Z"
    },
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "84aca8ea777720d8",
    "outputId": "b6ec0860-3546-411f-9dcc-bed57e4e54a3"
   },
   "id": "84aca8ea777720d8",
   "execution_count": 22,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#初值\n",
    "#loss_ic(self, x_ic, U_ic):\n",
    "num_i_train = 201\n",
    "id_ic = np.random.choice(num_x, num_i_train, replace=False)\n",
    "irhoL, iuL, ipL = 10.0, 1.0, 0.2\n",
    "irhoR, iuR, ipR = 0.1, 1.0, 0.2\n",
    "rho_u_p_LR = [irhoL, iuL, ipL, irhoR, iuR, ipR]\n",
    "x_ic = x_grid[id_ic, 0][:, None]\n",
    "t_ic = t_grid[id_ic, 0][:, None]\n",
    "x_ic_train = np.hstack((t_ic, x_ic))\n",
    "U_ic = IC(x_ic_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "#tensor\n",
    "x_ic_train = torch.tensor(x_ic_train, dtype=dtype, requires_grad=True, device=device)\n",
    "U_ic = torch.tensor(U_ic, dtype=dtype, device=device)\n",
    "#loss_xic(self, t_ic, x_ic):\n",
    "t_ic_train = t_grid[id_ic, 0][:, None]\n",
    "x_ic = t_ic_train*1\n",
    "#tensor\n",
    "t_ic_train = torch.tensor(t_ic_train, dtype=dtype, device=device, requires_grad=True)\n",
    "x_ic = torch.tensor(x_ic, dtype=dtype, device=device)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:46.980811Z",
     "iopub.execute_input": "2024-09-24T08:23:46.981641Z",
     "iopub.status.idle": "2024-09-24T08:23:47.183851Z",
     "shell.execute_reply.started": "2024-09-24T08:23:46.981604Z",
     "shell.execute_reply": "2024-09-24T08:23:47.183052Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.040415Z",
     "start_time": "2024-11-13T03:05:15.025925Z"
    },
    "id": "bb1f3322bfae61ab"
   },
   "id": "bb1f3322bfae61ab",
   "execution_count": 23,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.056071Z",
     "start_time": "2024-11-13T03:05:15.042511Z"
    },
    "id": "986db6b6dffa053e"
   },
   "cell_type": "code",
   "source": [
    "#bc\n",
    "# 边界设置,all time and x=lb,x=ub\n",
    "x_lbc = x_grid[0, :][:, None]\n",
    "x_rbc = x_grid[-1, :][:, None]\n",
    "t_bc = t_grid[0, :][:, None]\n",
    "x_lbc_train = np.hstack([t_bc, x_lbc])\n",
    "x_rbc_train = np.hstack([t_bc, x_rbc])\n",
    "U_lbc = IC(x_lbc_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "U_rbc = IC(x_rbc_train, irhoL, iuL, ipL, irhoR, iuR, ipR)\n",
    "x_bc_train = np.vstack([x_lbc_train, x_rbc_train])\n",
    "U_bc = np.vstack([U_lbc, U_rbc])\n",
    "#\n",
    "x_bc_train = torch.tensor(x_lbc_train, dtype=dtype, requires_grad=True, device=device)\n",
    "U_bc = torch.tensor(U_rbc, dtype=dtype, device=device)"
   ],
   "id": "986db6b6dffa053e",
   "execution_count": 24,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#残差\n",
    "#loss_pde(self, x_int):\n",
    "num_f_train = 5000\n",
    "id_f = np.random.choice(num_x * num_t, num_f_train, replace=False)\n",
    "x_int = X[:, 0][id_f, None]\n",
    "t_int = T[:, 0][id_f, None]\n",
    "x_int_train = np.hstack((t_int, x_int))\n",
    "#rh\n",
    "x_int_l_train, x_int_r_train = x_left_right(x_int_train, dx)\n",
    "#tensor\n",
    "x_int_train = torch.tensor(x_int_train, dtype=dtype, requires_grad=True, device=device)\n",
    "x_int_l_train = torch.tensor(x_int_l_train, dtype=dtype, requires_grad=True, device=device)\n",
    "x_int_r_train = torch.tensor(x_int_r_train, dtype=dtype, requires_grad=True, device=device)\n"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:47.599759Z",
     "iopub.execute_input": "2024-09-24T08:23:47.600812Z",
     "iopub.status.idle": "2024-09-24T08:23:47.609794Z",
     "shell.execute_reply.started": "2024-09-24T08:23:47.600765Z",
     "shell.execute_reply": "2024-09-24T08:23:47.608821Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.071456Z",
     "start_time": "2024-11-13T03:05:15.057074Z"
    },
    "id": "b5ad1a076b1101fa"
   },
   "id": "b5ad1a076b1101fa",
   "execution_count": 25,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "#实例化神经网络\n",
    "layers = Dnnlayers(2, 3, 5, 30)\n",
    "# Pinns(layers)\n",
    "model = Pinns(layers)\n",
    "loss_total_history = []\n",
    "loss_pde_history = []\n",
    "loss_ibc_history = []\n",
    "loss_rh_history = []\n",
    "lr = 0.001\n",
    "optimizer = torch.optim.Adam((model.dnn.parameters()), lr=lr)"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:49.532299Z",
     "iopub.execute_input": "2024-09-24T08:23:49.533197Z",
     "iopub.status.idle": "2024-09-24T08:23:50.456244Z",
     "shell.execute_reply.started": "2024-09-24T08:23:49.533152Z",
     "shell.execute_reply": "2024-09-24T08:23:50.455427Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.117853Z",
     "start_time": "2024-11-13T03:05:15.072458Z"
    },
    "id": "70089af1370df782"
   },
   "id": "70089af1370df782",
   "execution_count": 26,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-13T03:05:15.133063Z",
     "start_time": "2024-11-13T03:05:15.119856Z"
    },
    "id": "e103387bd2d6aa4e"
   },
   "cell_type": "code",
   "source": [
    "#training\n",
    "def train(epoch):\n",
    "    model.dnn.train()\n",
    "    def closure():\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        loss_pde = model.loss_pde(x_int_train)\n",
    "        loss_ic = model.loss_ic(x_ic_train, U_ic)\n",
    "        loss_bc = model.loss_bc(x_bc_train, U_bc)\n",
    "        loss_rh = model.loss_rh(x_int_train, x_int_l_train, x_int_r_train)\n",
    "\n",
    "        loss = loss_pde + 10*(loss_ic + loss_rh) # + loss_bc\n",
    "        if epoch%1==0:\n",
    "            print(f'epoch:{epoch}, loss:{loss:.8f},loss_pde:{loss_pde:.8f}, loss_IC:{loss_ic:.8f}，loss_BC:{+loss_bc:.8f} ,loss_rh:{loss_rh:.8f}')\n",
    "\n",
    "\n",
    "        loss_total_history.append(to_numpy(loss))\n",
    "        loss_pde_history.append(to_numpy(loss_pde))\n",
    "        loss_ibc_history.append(to_numpy(loss_ic))\n",
    "        loss_rh_history.append(to_numpy(loss_rh))\n",
    "\n",
    "        loss.backward()\n",
    "        return loss\n",
    "    loss = optimizer.step(closure)"
   ],
   "id": "e103387bd2d6aa4e",
   "execution_count": 27,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "epochs1 = 15000\n",
    "#for training\n",
    "tic = time.time()\n",
    "for epoch in range(1, epochs1+1):\n",
    "    train(epoch)\n",
    "toc = time.time()"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T08:23:52.150100Z",
     "iopub.execute_input": "2024-09-24T08:23:52.151126Z",
     "iopub.status.idle": "2024-09-24T09:31:41.292026Z",
     "shell.execute_reply.started": "2024-09-24T08:23:52.151082Z",
     "shell.execute_reply": "2024-09-24T09:31:41.291134Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:20:04.414442Z",
     "start_time": "2024-11-13T03:05:15.134193Z"
    },
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "59e6753e762b2257",
    "outputId": "7daa6d6d-ffed-4131-e36a-3cbaf6b11137"
   },
   "id": "59e6753e762b2257",
   "execution_count": 28,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "optimizer = torch.optim.LBFGS((model.dnn.parameters()),lr=1.0,max_iter=10)\n",
    "epochs2 =1500\n",
    "tic = time.time()\n",
    "for epoch in range(15001, 15000+epochs2+1):\n",
    "    train(epoch)\n",
    "toc = time.time()"
   ],
   "metadata": {
    "execution": {
     "iopub.status.busy": "2024-09-24T09:31:41.293728Z",
     "iopub.execute_input": "2024-09-24T09:31:41.294028Z",
     "iopub.status.idle": "2024-09-24T09:38:28.107433Z",
     "shell.execute_reply.started": "2024-09-24T09:31:41.293996Z",
     "shell.execute_reply": "2024-09-24T09:38:28.106494Z"
    },
    "trusted": true,
    "ExecuteTime": {
     "end_time": "2024-11-13T03:22:06.064738Z",
     "start_time": "2024-11-13T03:20:04.415445Z"
    },
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "cdd6739e70b2ed04",
    "outputId": "31180b13-1c61-404f-cf74-4f31a6142da1"
   },
   "id": "cdd6739e70b2ed04",
   "execution_count": 32,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-13T03:22:06.632203Z",
     "start_time": "2024-11-13T03:22:06.065741Z"
    },
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 452
    },
    "id": "5d53fe4a435204b4",
    "outputId": "debd01f4-ec04-45f5-daa1-0a19d3f7b162"
   },
   "cell_type": "code",
   "source": [
    "eporchx = np.arange(1, len(loss_total_history)+1)\n",
    "plt.plot(eporchx, loss_total_history, label='loss_total_history')\n",
    "plt.plot(eporchx, loss_pde_history, label='loss_pde_history')\n",
    "plt.plot(eporchx, loss_ibc_history, label='loss_ibc_history')\n",
    "plt.plot(eporchx, loss_rh_history, label='loss_rh_history')\n",
    "plt.yscale('log')\n",
    "plt.title('loss with eproch')\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "5d53fe4a435204b4",
   "execution_count": 30,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "data = scipy.io.loadmat(r'1Dexact_Case1.mat')  # Import Solution data\n",
    "x_test = data['x'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "t_test = data['t'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "exact_rho = np.real(data['rho'])  # rho(t,x,y)\n",
    "exact_p = np.real(data['p'])\n",
    "exact_u = np.real(data['u'])  # Exact total\n",
    "t_grid_test, x_grid_test = np.meshgrid(t_test, x_test)\n",
    "x_int_test = np.hstack((t_grid_test.flatten()[:, None], x_grid_test.flatten()[:, None]))\n",
    "x_int_test = torch.tensor(x_int_test, dtype=dtype, requires_grad=True, device=device)\n",
    "\n",
    "rho_predict, u_predict, p_predict = model.forword(x_int_test)\n",
    "t_grid_test = to_numpy(x_int_test[:, 0:1]).reshape(num_x, num_t).T\n",
    "x_grid_test = to_numpy(x_int_test[:, 1:2]).reshape(num_x, num_t).T\n",
    "rho_predict = to_numpy(rho_predict).reshape(num_x, num_t).T\n",
    "u_predict = to_numpy(u_predict).reshape(num_x, num_t).T\n",
    "p_predict = to_numpy(p_predict).reshape(num_x, num_t).T\n",
    "\n",
    "error_rho = abs(rho_predict - exact_rho)\n",
    "error_u = abs(u_predict - exact_u)\n",
    "error_p = abs(p_predict - exact_p)\n",
    "#\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, exact_rho, levels=np.linspace(exact_rho.min(), exact_rho.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_exact')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, rho_predict, levels=np.linspace(rho_predict.min(), rho_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_predict')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_predict, levels=np.linspace(u_predict.min(), u_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_predict')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_predict, levels=np.linspace(p_predict.min(), p_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_predict')\n",
    "plt.show()\n",
    "#\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_rho, levels=np.linspace(error_rho.min(), error_rho.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_rho')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_u, levels=np.linspace(error_u.min(), error_u.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_u')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_p, levels=np.linspace(error_p.min(), error_p.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_p')\n",
    "plt.show()\n",
    "\n",
    "print('L2 for error_rho is :', np.sqrt(((error_rho)**2).mean()))\n",
    "print('L2 for error_u is :', np.sqrt(((error_u)**2).mean()))\n",
    "print('L2 for error_p is :', np.sqrt(((error_p)**2).mean()))\n",
    "\n",
    "\n",
    "scipy.io.savemat(f'PINN_1Dpredict_Case1.mat', {'x': x_test, 't': t_test, 'rho': rho_predict, 'u': u_predict, 'p': p_predict})\n"
   ],
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-13T03:22:12.917557Z",
     "start_time": "2024-11-13T03:22:06.633220Z"
    },
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "f56e33350875c610",
    "outputId": "a247562e-5bc5-413a-fd56-7475e273a3a2"
   },
   "id": "f56e33350875c610",
   "execution_count": 31,
   "outputs": []
  },
  {
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "c62da339a82b0456",
    "outputId": "c99b3ee0-c085-4d94-be4a-2de64dff9f57",
    "ExecuteTime": {
     "end_time": "2024-12-03T13:31:25.078366Z",
     "start_time": "2024-12-03T13:31:17.173855Z"
    }
   },
   "cell_type": "code",
   "source": [
    "data = scipy.io.loadmat(r'1Dexact_Case1.mat')  # Import Solution data\n",
    "x_test = data['x'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "t_test = data['t'].flatten()[:, None]  # Partitioned spatial coordinates\n",
    "rho_exact = np.real(data['rho'])  # rho(t,x,y)\n",
    "p_exact = np.real(data['p'])\n",
    "u_exact = np.real(data['u'])  # Exact total\n",
    "data = scipy.io.loadmat(r'PINN_1Dpredict_Case1.mat')  # Import Solution data\n",
    "rho_predict = np.real(data['rho'])\n",
    "p_predict = np.real(data['p'])\n",
    "u_predict = np.real(data['u'])\n",
    "\n",
    "t_grid_test, x_grid_test = np.meshgrid(t_test, x_test)\n",
    "x_int_test = np.hstack((t_grid_test.flatten()[:, None], x_grid_test.flatten()[:, None]))\n",
    "x_int_test = torch.tensor(x_int_test, dtype=dtype, requires_grad=True, device=device)\n",
    "\n",
    "\n",
    "\n",
    "error_rho = abs(rho_predict - rho_exact)\n",
    "error_u = abs(u_predict - u_exact)\n",
    "error_p = abs(p_predict - p_exact)\n",
    "#\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, rho_exact.T, levels=np.linspace(rho_exact.min(), rho_exact.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_exact')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_exact.T, levels=np.linspace(u_exact.min(), u_exact.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_exact')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_exact.T, levels=np.linspace(p_exact.min(), p_exact.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_exact')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, rho_predict.T, levels=np.linspace(rho_predict.min(), rho_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('rho_predict')\n",
    "plt.savefig('PINNs1Drho_predict_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, u_predict.T, levels=np.linspace(u_predict.min(), u_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('u_predict')\n",
    "plt.savefig('PINNs1Du_predict_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, p_predict.T, levels=np.linspace(p_predict.min(), p_predict.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('p_predict')\n",
    "plt.savefig('PINNs1Dp_predict_Case1.png')\n",
    "plt.show()\n",
    "#\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_rho.T, levels=np.linspace(error_rho.min(), error_rho.max()+0.05, 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_rho')\n",
    "plt.savefig('PINNs1Derror_rho_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_u.T, levels=np.linspace(error_u.min(), error_u.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_u')\n",
    "plt.savefig('PINNs1Derror_u_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.contourf(x_grid_test, t_grid_test, error_p.T, levels=np.linspace(error_p.min(), error_p.max(), 1000), cmap='jet')\n",
    "plt.colorbar()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.title('error_p')\n",
    "plt.savefig('PINNs1Derror_p_Case1.png')\n",
    "plt.show()\n",
    "\n",
    "print('L2 for error_rho is :', np.sqrt(((error_rho)**2).mean()))\n",
    "print('L2 for error_u is :', np.sqrt(((error_u)**2).mean()))\n",
    "print('L2 for error_p is :', np.sqrt(((error_p)**2).mean()))\n",
    "'''\n",
    "L2 for error_rho is : 0.03384217426713135\n",
    "L2 for error_u is : 0.00032142902352134056\n",
    "L2 for error_p is : 0.00016079423958703805\n",
    "'''"
   ],
   "id": "c62da339a82b0456",
   "execution_count": 4,
   "outputs": []
  },
  {
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 843
    },
    "id": "25c168ff46c1fdb8",
    "outputId": "0e2723c8-cf6b-4d50-d320-5c16e8439b29",
    "ExecuteTime": {
     "end_time": "2024-12-03T13:31:26.593113Z",
     "start_time": "2024-12-03T13:31:26.397881Z"
    }
   },
   "cell_type": "code",
   "source": [
    "plt.figure()\n",
    "plt.plot(x_test, rho_exact[-1,:], 'b', label='end_rho_exact')\n",
    "plt.scatter(x_test, rho_predict[-1,:], color='r', s=1, label='end_rho_predict')\n",
    "plt.legend()\n",
    "#plt.savefig('1Drho_end_Case1.png')\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.plot(x_test[10:150], rho_exact[-1,10:150], 'b', label='end_rho_exact')\n",
    "plt.scatter(x_test[10:150], rho_predict[-1,10:150], color='r', s=1, label='end_rho_predict')\n",
    "plt.legend()\n",
    "#plt.savefig('1Drho_end_Case1.png')\n",
    "plt.show()"
   ],
   "id": "25c168ff46c1fdb8",
   "execution_count": 5,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Ko8Sdd_KpHY-"
   },
   "id": "Ko8Sdd_KpHY-",
   "execution_count": null,
   "outputs": []
  }
 ]
}
